
def LUd(a): #a is matrix
  #use with LUbacksub to solve or invert
  print 'a starts:',a
  n=len(a)
  indx=[0 for i in range(n)]
  i=imax=j=k=0
  big=dum=sm=temp=0.0
  
  vv=[0.0 for i in range(n)] #stores implicit scaling for each row
  d=1.0 #no row interchanges yet
  #get implicit scaling information
  for i in range(len(a)): #loop rows
    big=0.0 #stores maximum value found
    for j in range(len(a)): #loop columns
      if abs(a[i][j])>big:big=a[i][j] #store biggest of the two
      vv[i]=1.0/big #stores scaling factor
  #loop columns for Crout's method
  for j in range(n): #loop columns
    for i in range(j-1): #loop rows up to current column
      sm=a[i][j] #start sum off with entry
      for k in range(i-1): #loop rows until reach current row number
        sm-=a[i][k]*a[k][j] #subtract row,
        a[i][j]=sm
    big=0.0
    for i in range(j,n):
      sm=a[i][j]
      for k in range(j):
        sm-=a[i][k]*a[k][j]
      a[i][j]=sm
      dum=vv[i]*abs(sm)
      if dum>=big:
        big=dum
        imax=i
    if j!=imax: #need to interchange rows
      for k in range(n):
        dum=a[imax][k]
        a[imax][k]=a[j][k]
        a[j][k]=dum
      d=-d
      vv[imax]=vv[j]
    indx[j]=imax
    if a[j][j]==0.0:a[j][j]=10**-20 #small number
    if j!=n: #divide by pivot element
      dum=1.0/a[j][j]
      for i in range(j+1,n):
        a[i][j]*=dum
    print 'a ends:',a
    return a,indx,d

def LUbacksub(a,indx,b):
  ii=0
  sm=0.0
  n=len(a)
  print b

  for i in range(n):
    print 'i=',i
    ip=int(indx[i])
    sm=b[ip]
    b[ip]=b[i]
    if ii!=0:
      for j in range(ii,i-1):
        print j
        sm-=a[i][j]*b[j]
    elif sm!=0:ii=i
    b[i]=sm
  print 'first loop',b
  for i in range(n-1,0,-1):
    sm=b[i]
    for j in range(i+1,n):
      sm-=a[i][j]*b[j]
    b[i]=sm/a[i][i]
    print 'stored',b[i]
  return b


def doLUsolve(a,b):
  [a,indx,d]=LUd(a)
  x=LUbacksub(a,indx,b)
  return x
